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Thermal boundary conditions has played an increasingly important role in revealing the nature 
of short-range spin glasses and is likely to be relevant also for other disordered systems. Diffusion 
method initializing each replica with a random boundary condition at the infinite temperature 
using population annealing has been used in recent large-scale simulations. However, the efficiency 
of this method can be greatly suppressed because of temperature chaos. For example, most samples 
have some boundary conditions that are completely eliminated from the population in the process of 
annealing at low temperatures. In this work, I study a weighted average method to solve this problem 
by simulating each boundary conditions separately and collect data using weighted averages. The 
efficiency of the two methods are studied using both population annealing and parallel tempering, 
showing that the weighted average method is more efficient and accurate. 


I. INTRODUCTION 

Thermal boundary conditions (TBC) includes the set 
of all combinations of periodic/antiperiodic boundary 
conditions in each spatial dimension m- For exam¬ 
ple in three dimensions (d = 3), there are 2“^ = 8 possible 
choices. Each boundary condition i occurs in the ther¬ 
mal ensemble with weight pi depends on its free energy 
Fi as Pi = exp[—/3(Fi — F)], where F is the total free en¬ 
ergy of the system in TBC and (3 is the inverse temper¬ 
ature. TBC was initially motivated to reduce domain- 
wall effects in spin glasses [T] and was later shown to 
be useful in studying temperature chaos and bond chaos 
via boundary condition crossings, i.e, the weights {pi} 
change chaotically as a function of temperature or cou¬ 
plings la Previously, thermal boundary conditions 
was used with exact algorithms for finding ground states 
of two-dimensional spin glasses [U [S] (referred to as “ex¬ 
tended” boundary conditions). A restricted version of 
thermal boundary conditions using periodic and anti- 
periodic boundary conditions in a single direction (a sub¬ 
set of TBC) was also used in Refs. [BHS]- 

Simulating the full TBC ensemble in Refs. m was 
done using population annealing (PA) [TDHT5] . See 
Ref. [13] for a recent discussion of the algorithm. Popu¬ 
lation annealing works with a population of replicas and 
aims to maintain thermal equilibrium while lowering the 
temperature. When temperature is decreased, the popu¬ 
lation is resampled according to the Boltzmann weights 
of the replicas, followed by regular Monte Carlo sweeps to 
all replicas in the population. In this work, the Metropo¬ 
lis algorithm is used. Simulating TBC using the diffusion 
method with PA is very simple. One can initialize each 
replica with a random boundary condition at /3 = 0 and 
then boundary conditions are resampled along with the 
replicas. The word “diffusion” would become apparent 
when one implements the method in parallel tempering. 
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as discussed in Sec. Eel 

It was noticed in Refs. [HE] that the efficiency of the 
diffusion method can be greatly suppressed by temper¬ 
ature chaos. Some boundary conditions can be totally 
removed from the population. This is not satisfactory 
as there is no mechanism to recover these lost bound¬ 
ary conditions once they are eliminated from the pop¬ 
ulation. Furthermore, temperature chaos predicts that 
these boundary conditions could become important again 
at lower temperatures. Therefore, it is worth to study a 
new method that does not have this problem. The most 
natural way is to simulate each boundary condition sep¬ 
arately and combine the observables using weighted av¬ 
erages with free energy. How to weight different kinds of 
observables is discussed in Sec. umi Since it is also in¬ 
teresting to study the performance of parallel tempering 
(PT) in simulating the full set of thermal boundary con¬ 
ditions, the two methods are therefore studied in both 
PA and PT. 

The paper is organized as follows: Section|ll]introduces 
the Edwards-Anderson model, measured quantities, the 
diffusion method and the weighted average method. Nu¬ 
merical results are shown in Sec. followed by conclud¬ 
ing remarks and future challenges in Sec. |IV| 

II. MODELS, OBSERVABLES AND METHODS 
A. The Edwards-Anderson model 

The Edwards-Anderson (EA) Hamiltonian is defined 
as 

H = — JijSiSj, (1) 

iu) 

where {Si = ±1} are the spin degrees of freedom de¬ 
fined on a three-dimensional cubic lattice. The sum over 
(ij) means sum over all nearest neighbours. Jij is the 
coupling between spins Si and Sj and all couplings are 
independently chosen from the standard Gaussian distri¬ 
bution with mean 0 and variance 1. Thermal boundary 
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conditions is applied to all instances. 


B. Observables 

There are three classes of observables that need to be 
treated differently using weighted averages. The first 
class are observables that are functions of a single replica 
like the energy E or the magnetization m. The second 
class are the thermodynamic observables of the free en¬ 
ergy F and the entropy S. There is a third class of ob¬ 
servables in spin glasses due to the nature of the symme¬ 
try breaking which are functions of two replicas like the 
spin overlap q defined as 

(2) 

i 

where micro-states a, b are chosen independently from 
the Boltzmann distribution. Note that a, b can be chosen 
from the same boundary conditions or different bound¬ 
ary conditions. Another quantity in this class is the link 
overlap which needs some care in the definition due to 
the change of boundary conditions, the link overlap qi is 
defined as 

^ E 5“sign( J“ )^“5^ign(4)*5l, (3) 

(b> 

where the sign function is ±1 depends on whether the 
argument is positive or negative respectively. The defini¬ 
tion has nothing to do with the weighted average method, 
but this is a more useful definition in the TBC setting 
for studying the nature of spin glasses. In this way, the 
difference of the link overlap of two different boundary 
conditions arises only from the domain walls, not from a 
mixture of domain walls and boundary condition differ¬ 
ences. In the following, I will focus on the spin overlap 
function and will refer to it as the overlap function where 
no confusion arises. Note also that there is a 8 x 8 over¬ 
lap function matrix in the TBC ensemble. An important 
statistic /, which quantifies the weight of an overlap dis¬ 
tribution near g = 0, is defined as 

1-0.2 

1= / Piq)dq, (4) 

J-0.2 

where P{q) is an overlap distribution function. 


C. The diffusion method 

We have already discussed the diffusion method in pop¬ 
ulation annealing, which initializes each replica with a 
random boundary condition at /3 = 0. For completeness, 
I also introduce and study the diffusion method in paral¬ 
lel tempering because parallel tempering is widely used 
in spin-glass simulations. Following the same strategy, 
one can simulate the TBC using parallel tempering by 


generating random states with random boundary condi¬ 
tions at /3 = 0. Then the boundary conditions diffuse 
along with the replicas under the swap moves of parallel 
tempering, hence the name the diffusion method. The 
convenience of working with /3 = 0 is because propos¬ 
ing a boundary condition change at finite temperatures 
can be costly as many bonds are affected when boundary 
condition changes. Furthermore, this is also essential to 
measure the absolute free energy. 

The implementation of this method is simple and de¬ 
tailed balance is satisfied. However, the efficiency of the 
method still needs to be studied. Since in this method, 
PA and PT work in a similar way, one may expect that 
PT should suffer from temperature chaos too. This turns 
out to be the case as discussed in Secs. IIII Al and IIIIBI 
The effect of boundary condition eliminations in PA is 
replaced by one or more diffusion bottlenecks in PT. In 
the next, I discuss the weighted average method, which 
does not have this problem. 


D. The weighted average method 


It was shown that the absolute free energy can be mea¬ 
sured very accurately using the free energy perturbation 
method in both population annealing and parallel tem¬ 
pering miiin]. It can be shown from statistical mechan¬ 
ics that the average energy, entropy, free energy and the 
spin overlap distribution should be averaged as 


W’ 

II 

(5) 

s = y^^s^pi -y^^Pi log Pi, 

i i 

(6) 

p = y2 ^ E ^ogpi, 

i i 

(7) 

p(.<i) = y2Piji<i)p^pj’ 

(8) 




where pi 


e-PP 
^ e-PE ’ 


Ei , Si and Fi are the energy, en¬ 


tropy and free energy of boundary condition i and (q) 
is the overlap distribution between the boundary condi¬ 
tions i and j. 

It is worth noting that the weighted average method 
generates a lot more information than the diffusion 
method. For example, the overlap matrix is briefly dis¬ 
cussed in Sec. HYl For now, we discuss briefly of the 
implementation of the weighted average method in PA 
and PT. In the weighted average method, each boundary 
condition is simulated independently. It is only when 
computing the spin overlaps that one needs communi¬ 
cations between replicas with different boundary condi¬ 
tions. Therefore, in this way, PT can be simulated using 
parallel computing too. Since it is usually a practice to 
simulated two independent sets of replicas of each bound- 
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ary condition, one can use 8 or 16 threads in the simula¬ 
tion. This is not doable in the diffusion method of PT. 
On the other hand, PA is intrinsically parallel and is also 
much more flexible with the number of threads. In my 
simulations, I have used OpenMP parallel computing for 
both PA and PT. For the equilibration measure of PA, 
one can either use the weighted average or the minimum 
of the entropy of families m of all boundary conditions. 


III. RESULTS 


In this section, I investigate the performance of the 
diffusion method and the weighted average method. I 
first compare PA and PT for the diffusion method in 
Sec. IIII Al The conclusion is that the diffusion method, 


while works but can be inefficient and the limitations of 
the method are discussed in Sec. |IIIB] Finally, I study the 
performance of the weighted average method, compare 
PA and PT and also with the diffusion method, showing 
that the weighted average method works better than the 
diffusion method in Sec. IIII Cl 


TABLE I: The simulation parameters of PA using the 
diffusion method (DF) for different system sizes L with ther¬ 
mal boundary conditions. R represents the total number of 
replicas, l//9max is the lowest temperature simulated, Nt the 
number of temperatures in the annealing schedule evenly dis¬ 
tributed in p, Ns the number of sweeps per replica per tem¬ 
perature and M the number of samples. The simulation pa¬ 
rameters using the weighted average method (WA) are also 
shown in the table, where R represents the population size of 
each boundary condition. 


Method 

L 

R 

1 / /3max 

Nt 

Ns 

M 

DF 

4 

510-* 

0.2 

101 

10 

101 

DF 

6 

210® 

0.2 

101 

10 

101 

DF 

8 

510® 

0.2 

201 

10 

101 

DF 

10 

210® 

0.2 

301 

10 

101 

WA 

4 

510® 

0.2 

101 

10 

101 

WA 

6 

210'* 

0.2 

101 

10 

101 

WA 

8 

510-* 

0.2 

201 

10 

101 

WA 

10 

210® 

0.2 

301 

10 

101 

WA 

12 

210® 

1/3 

301 

10 

101 


A. PA and PT: the diffusion method 

In this section, I compare the efficiency of the diffusion 
method of PA and PT in simulating thermal boundary 
conditions. The large-scale population annealing data is 
taken from a recent work on temperature chaos [5] . In the 
implementation of parallel tempering, I used Njs temper¬ 
atures evenly distributed in /3 for the high temperature 
part and temperatures evenly distributed in T for the 
low temperature part. No Monte Carlo sweeps were done 
at /3 = 0, it is sufficient to generate random states and 
boundary conditions at the infinite temperature. The 
amount of work in this work is counted in Monte Carlo 
sweeps and one Monte Carlo sweep is a sequential up¬ 
date of all the spins for a replica at a temperature once. 
The work is distributed evenly between thermalization 
and data collection. The simulation parameters of PA 
and PT are summarized in Tables and [ll] respectively. 
The main conclusion of this section is that PA and PT 
has similar efficiency in the diffusion method. The limi¬ 
tations of the method are discussed in Sec. IIII Bl 

1. Detailed comparison of a single hard sample 

In this section, I make a detailed comparison for the 
hardest instance of L = 8 chosen from a total of 2001 
instances in Ref. [2]. I will focus on the systematic error 
of —(3F, the q distribution at a low temperature /3 = 5 as 
well as the evolution of {pi} as a function of /3. One can 
see from the evolution of {pi} that this sample is rather 
chaotic with multiply boundary condition crossings. 

The convergence of the dimensionless quantity —j3F 
at a low temperature /3 = 5 is studied as a function 


TABLE II: The simulation parameters of PT using the dif¬ 
fusion method (DF) for different system sizes L with thermal 
boundary conditions. Ng represents the number of temper¬ 
atures in the high temperature part uniformly distributed in 
P, Nt the number of temperatures in the low temperature 
part uniformly distributed in T, Tmin the lowest temperature 
simulated, Tmax the highest temperature simulated in the low 
temperature part. Ns the number of sweeps per temperature 
and M the number of samples. The work is evenly distributed 
in two independent sets of replicas. The simulation parame¬ 
ters using the weighted average method (WA) are also shown 
in the table, where Ns represents the number of sweeps of 
each boundary condition. 


Method 

L 

Ng 

Nt 

Tmin 

Tmax 

Ns 

M 

DF 

4 

5 

10 

0.2 

2.0 

210® 

101 

DF 

6 

5 

20 

0.2 

2.0 

410® 

101 

DF 

8 

10 

30 

0.2 

2.0 

1.210^ 

101 

DF 

10 

10 

40 

0.2 

2.0 

610^ 

101 

WA 

4 

5 

10 

0.2 

2.0 

210® 

101 

WA 

6 

5 

20 

0.2 

2.0 

410® 

101 

WA 

8 

10 

30 

0.2 

2.0 

1.210® 

101 

WA 

10 

10 

40 

0.2 

2.0 

610® 

101 

WA 

12 

20 

40 

1/3 

2.0 

510® 

101 


of the amount of work W. Even though two indepen¬ 
dent sets of replicas are simulated in parallel tempering 
to compute the overlap distribution and data were col¬ 
lected from both sets, when studying the convergence of 
the free energy, this factor of 2 is not counted since this 
only improves statistical errors, not systematic errors. 
Indeed, measuring free energy does not require two sets 
of independent runs. One can say that it is an advantage 
of population annealing that can effectively compute the 
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overlap distribution from one run due to the ensemble 
nature of the algorithm. The same parameters in Ta¬ 
bles U and |TT1 are used for this instance and the amount of 
work is varied by changing the population size in PA and 
the number of sweeps in PT. The convergence of —j3F is 
shown in the top panel of Fig. [l] The standard is taken 
from a weighted average m of the largest data set of PA, 
and the errorbars in the plot are estimated from multiple 
independent runs. The overlap distribution and the evo¬ 
lution of {pi \ as a function of /3 are shown in Fig.[^ The 
results are randomly chosen from the largest runs of PA 
and PT in Fig. We see that PA and PT have similar 
performance in the diffusion method. 



logi„VV 

FIG. 1: (Color online) Linear-log plot of the systematic error 
A(—/JF) of PA and PT as a function of the amount of work 
W for a chaotic instance of L = 8 at T = 0.2. Top panel: the 
diffusion method (DF). Bottom panel: the weighted average 
method (WA). The errorbars are the standard deviation of 
the —PF distribution computed from multiple independent 
runs, not the errorbar of the sample mean of —fiF. Note 
that even though the weighted average simulates all the 8 
boundary conditions independently, it still converges faster 
than the diffusion method and is therefore more efficient. 


2. A large-scale comparison 

Now we look at a large-scale comparison of PA and PT 
using the diffusion method. Two important quantities for 
spin glasses are compared: the free energy per spin / and 
the statistic I of the overlap distribution. The results for 
/ and / are shown in Figs. and respectively, again 
showing the similar performance of PA and PT. 


B. Limitations of the diffusion method 

We have shown in the last section that one can use 
the diffusion method to simulate thermal boundary con¬ 
ditions. This method is simple to implement and yields 
fair result. However, the method has its own limitations 
due to the diffusion nature of the method, in particular, 
the problem of the boundary condition dying out. 



0 1 2 3 4 5 

P 

FIG. 2: (Color online) The comparison of the evolution of {pi} 
as a function of P of PA and PT using the diffusion method. 
The solid lines are for population annealing and the points 
are for parallel tempering. Each color represents a boundary 
condition. The inset is the comparison of the overlap distri¬ 
bution at d = 5. The results are randomly chosen from the 
largest runs of Fig. 
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/pa 

FIG. 3: (Color online) Comparison of the free energy per spin 
/ at (3 = 5 of population annealing and parallel tempering 
using the diffusion method. 


The problem of boundary condition dying out of the 
diffusion method is actually a very severe one. In the 
state-of-the-art simulation of Ref. [2], where most sam¬ 
ples if not all are well equilibrated with decent simulation 
parameters, it is striking that most of the samples suffer 
from this problem as shown in Table and Fig. One 
can see that about 77% of instances have some boundary 
conditions eliminated at T =1/3 and the number quickly 
grows to about 98% at T = 0.2. The distribution of the 
number of boundary conditions that are eliminated no 
also quickly shifts from a bias towards nu = 0 to = 7. 

This problem is a resolution problem as the resolution 
of the weight of a boundary condition is 1/R. A bound- 
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FIG. 4: (Color online) Comparison of the statistic I at P = 
5 of population annealing and parallel tempering using the 
diffusion method. 
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FIG. 5: (Color online) Two typical distribution of the num¬ 
ber of boundary conditions that die out n_D of L = 10 at 
T = 1/3 and 0.2 respectively in the diffusion method. Note 
that as temperature is lowered a little in the low-temperature 
phase, a dramatic change occurs in the distribution as many 
boundary conditions are eliminated in the annealing process. 
In contrast, no boundary condition dies out in the weighted 
average method. 


TABLE III: The fraction of instances with boundary con¬ 
dition dying out fj in the large-scale simulation of Ref. 
with M — 2001 samples of each system size L. In contrast, no 
boundary condition dies out in the weighted average method. 


L 

4 

6 

8 

10 

10 

12 

T 

0.2 

0.2 

0.2 

0.2 

1/3 

1/3 

fj 

0.9835 

0.9840 

0.9820 

0.9720 

0.7701 

0.7731 


ary condition i is likely to be eliminated when pi becomes 
close to 1/i?: pi ~ 1/i?. When pi < 1/R, it is not even 
expected to be part of the ensemble in most independent 
runs. Unluckily, there is no mechanism in the diffusion 
method to reintroduce a boundary condition once it is 
eliminated. To “confirm” that the simulation is fine, one 
can only try a much larger i?, like a factor of 10 as used 
in Ref. [2], and see if the results remain the same. But 
this is not satisfactory as the runs can be consistently 
wrong and this is also not efficient. This is especially the 
case for large system sizes as temperature chaos asserts 
that the dominate boundary conditions switches more 
frequently as system size increases. The problem of the 
total absence of some boundary conditions also occurs in 
the parallel tempering implementation. Bottlenecks in 
the diffusion channel prevents some boundary conditions 
to reach the lowest temperature. This is not too surpris¬ 
ing considering the similar performance of PA and PT. 
Therefore, this is a limitation of the diffusion method 
itself, independent of which algorithm one uses to imple¬ 
ment it. In the next section, I study the weighted average 
method, its validity as well as its efficiency. 


C. PA and PT: the weighted average method 

In this section, I present a detailed study of the same 
hard instance as well as a large-scale simulation using 
the weighted average method with PA and PT. At first 
sight, this might be less efficient as one has to simulate 
all the 8 boundary conditions independently. However, 
as we shall see in this section that this is not the case. 
For the hard instance studied, it is even more efficient 
than the diffusion method. In addition, the problem of 
boundary conditions elimination is completely removed 
in the weighted average method. 


1. Detailed comparison of the single hard sample 

The convergence of —/fF is shown in the bottom panel 
of Fig. [2 Compare with the diffusion method in the top 
panel, one can see that it is more efficient. This is in¬ 
teresting and therefore the weighted average method is 
both more efficient and more accurate. PA and PT are 
still comparable in performance. One may have noticed 
that PT converges slightly faster, but there are more fac¬ 
tors to consider here. One is that PT often requires two 
independent sets of replicas, which are not counted here, 
to compute the spin overlaps. Take this into account, 
PA and PT are similar in performance. Another factor is 
that PT also has slightly more overhead in parallel com¬ 
puting due to the frequent breaks to measure the overlap 
functions. However, one can choose to run the code in 
sequential like in the diffusion method. Therefore, it is 
more fair to conclude that PA and PT are comparable in 
efficiency. The overlap distribution and the evolution of 
{pi} as a function of /3 are shown in Fig. The results 
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are randomly chosen from the largest runs of PA and PT 
in Fig. The two algorithms give very similar result. 



P 


FIG. 6: (Color online) The comparison of the evolution of 
{pi} as a function of P of PA and PT using the weighted 
average method. The solid lines are for population annealing 
and the points are for parallel tempering. The results are 
randomly chosen from the largest runs of Fig. The inset is 
the comparison of the overlap distribution at /3 = 5. Note that 
even with slightly less computational work, the agreement is 
better than that of Fig. [fusing the diffusion method, showing 
that the weighted average method is more accurate. 


2. A large-scale simulation 

Finally, we do a large-scale comparison of PA and PT. 
We have studied the EA model up to L = 10 at T = 0.2 
and L = 12 at T =1/3 (current state-of-the-art) using 
the weighted average method, demonstrating that the 
method can be used for an accurate large-scale study 
of spin glasses with thermal boundary conditions. The 
comparison of the free energy per spin / and the statistic 
I of the overlap distribution are shown in Figs. and 
respectively, again showing the similar performance of 
PA and PT. It also worth noting that the agreement in 
the measurements of I is better than that the diffusion 
method even though with slightly less amount of work, 
in line with the results of the single hard instance. 



/pa 

FIG. 7: (Color online) Comparison of the free energy per 
spin / of population annealing and parallel tempering using 
the weighted average method. P = 5 for L — 6,8 and 10 and 
/3 = 3 for L = 12. 
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FIG. 8: (Color online) Comparison of the statistic I of pop¬ 
ulation annealing and parallel tempering using the weighted 
average method. /? = 5 for L = 6,8 and 10 and /3 = 3 for 
L = 12. 


IV. CONCLUSIONS AND FUTURE 
CHALLENGES 

In this paper, I studied two methods to simulate ther¬ 
mal boundary conditions, the diffusion method and the 
weighted average method using both population anneal¬ 
ing and parallel tempering. I have illustrated the weak¬ 
nesses of the diffusion method, in particular, the elim¬ 
ination of boundary conditions. The weighted average 
method is studied in detail, showing that it is more ef¬ 
ficient than the diffusion method and solves the bound¬ 
ary condition elimination problem. Population annealing 


and parallel tempering on the other hand are similar in 
their performance. 

Furthermore, the weighted method method generates a 
lot more information than the diffusion method and has 
opened new ways in doing spin-glass research in ther¬ 
mal boundary conditions. Detailed information of each 
boundary condition is accessible, and the spin and link 
overlap distributions between any two boundary condi¬ 
tions are also readily available. For example, it is inter¬ 
esting to study the overlap distributions in the whole in- 
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stance as well as in a window between the same boundary 
conditions and two different boundary conditions. This 
can be very useful in studying the nature of domain walls, 
i.e. whether the low-lying excitations are space filling or 
have fractal surfaces. It is also interesting to study the 
structure of the spin overlap matrix. Except studying the 
size dependence of the traditional I in periodic bound¬ 
ary conditions and thermal boundary conditions, one can 
also study the weighted or minimum I of the diagonal 
part of the overlap matrix. In this way, one can remove 
effects of the overlaps between different boundary condi¬ 
tions, which are generally larger than the diagonal ones 
near the origin, and therefore may further reduce finite- 
size effects. Efforts along these directions are currently 


in progress and will be presented in future publications. 
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